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Interesting and recommended websites: 


Barry Kissane (Australia) has a great website which collects many promising links. I strongly 
recommend visiting: 

http://wwwstaff.murdoch.edu.au/-kissane/pd/internetmaths.htm 


You can download the Proceedings of CAME 2009. Thanks to Djordje Kadijevich from Serbia. 

http://www.megatrend.edu.rs/came_files/CAME%202009-Proceedings.pdf 


DUG-Member Rene Hugelshofer (Switzerland) informed about a new 65 pages pdf-paper on 
quadratics treated with TI-NspireCAS (in German): 


Quadratische Funktionen 
und Gleichungen mit CAS 

by Benno Frei, Rene Hugelshofer and Robert Marki 

x\ufgabe 4,14 


Vt)n oiner touron Marmorplattc i.st boi dor Boarbeitung 
ein btiick abgobrochen. Din Bruchlinie EF istoin Cjora- 
denstuLk. L'm don Schadon ini>glichst kloin /u halton, 
schnuidcn wir aus don boidon BruLhstiickon jew oils auf 
die skizzierte VVeise oino roditeekige Platte so aus^ dass 
die bumme der Flaeheninhalte mdglichst gross wild. 
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This paper - and many others can be downloaded from: 

http://www.ti-unterrichtsmaterialien.net 


DUG-Member Eberhard Lehmann published his 

Matrizenrechnung, Anwendungen Teil 2 

On 352 pages Eberhard demonstrates many applications of matrices reaching from mappings to 
Markov chains. Eberhard uses TI-Nspire, Voyage 200, DERIVE and his own software ANIMA- 
TO. More details can be found at his website 


http://home.snafu.de/mirza 

If you know about interesting websites and publications then please inform me for 
sharing the information in our community. Josef 
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Dear DUG Members, 

Welcome to DNL#75 - we con have a little celebration. We need "only" 25 more 
DNLs to reach the 1001! 

I'd like to welcome also a couple of new DUG-Members. It is great that our world 
wide community is still growing despite the fact that DERIVE is not available any 
longer. I was very happy receiving Danny Ross Lunsford's short note on a request 
which was published in the DERIVE-NEWS mailing list (see page 44). 

There are not so many contributions in this newsletter as usual. You will see that 
Michel's and Jose Luis' articles - both are presentations from ACA 09 - are very 
extended. It will need the next DNL, too, in order to present the missing papers. 

This DNL75-pdf-file has 2 MB because I included the ready made pdf-file of Jose 
Luis Golan's contribution on RANDOM SAMPLES. It would have been too difficult 
to rewrite it in WORD including all his great illustrations. 

In addition to these two papers I included another project for students provided by 
Roland Schroder. The Dog & Biker problem is a nice application of a recursive pro- 
cedure and gives the chance to demonstrate students the power of the ITERATES- 
function of DERIVE. 

At the other hand recursive procedures are very suitable for spreadsheet pro- 
grams. So I tried to transfer this problem to TT-NspireCAS and could implement 
some kind of animation. We - Roland Schroder and I - have still the problem to find 
the locus of the dog's jumps. It would be great if somebody could present the solu- 
tion. 

Tania Koller's valuable hint (User Forum) demonstrates that students of secondary 
schools are very familiar with hard- and software and they find solutions for occur- 
ing problems. Many thanks to Tania and her bright students. 

Best regards as ever 


Finally I'd like to remind you once more on our next conference which will be held in 
Andalusia, Spain. You can find how to submit and details about registration within 
the next few days. 

The Conference website is www.time2010.uma.es. 


Download all DA^L-Derive- and Tl-files from 

http : / / WWW . austromath . at/ dug/ 
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Contributions: 

Please send all contributions to the Editor. 
Non-English speakers are encouraged to 
write their contributions in English to rein- 
force the international touch of the DNL. It 
must be said, though, that non-English 
articles will be warmly welcomed nonethe- 
less. Your contributions will be edited but 
not assessed. By submitting articles the 
author gives his consent for reprinting it in 
the DNL. The more contributions you will 
send, the more lively and richer in contents 
the DERIVE & CAS-TI Newsletter will be. 


Next issue: December 2009 

Deadline 15 November 2009 

Preview: Contributions waiting to be published 

Some simulations of Random Experiments, J. Bohm, AUT, Lorenz Kopp, GER 

Wonderfui Worid of Pedai Curves, J. Bohm 

Toois for 3D-Problems, P. Liike-Rosendahl, GER 

Financial Mathematics 4, M. R. Phillips 

Hiil-Encription, J. Bohm 

Simuiating a Graphing Calcuiator in DERIVE, J. Bohm 
Henon, Mira, Gumowski & Co, J. Bohm 

Do you know this? Cabri & CAS on PC and Handheld, W. Wegscheider, AUT 

Steiner Point, P. Liike-Rosendahl, GER 

Overcoming Branch & Bound by Simuiation, J. Bohm, AUT 

Diophantine Poiynomiais, D. E. McDougaii, Canada 

Graphics Worid, Currency Change, P. Charland, CAN 

Cubics, Quartics - interesting features, T. Koller & J. Bohm 

Logos of Companies as an inspiration for Math Teaching 

Exciting Surfaces in the FAZ / Pierre Charland 's Graphics Gailery 

BooieanPlots.mth, P. Schofield, UK 

Old traditional examples for a CAS - what's new? J. Bohm, AUT 

Truth Tables on the Tl, M. R. Phillips 

Advanced Regression Routines for the TIs, M. R. Phillips 

Where oh Where is IT? (GPS with CAS), C. & P. Leinbach, USA 

Embroidery Patterns, H. Ludwig, GER 

Mandelbrot and Newton with DERIVE, Roman Hasek, CZ 

Snail-shells, Piotr Trebisz, GER 

A Conics-Explorer, J. Bohm, AUT 

Coding Theory for the Classroom?, J. Bohm, AUT 

Tutorials for the NSpireCAS, G. Herweyers, BEL 

Some Projects with Students, R. Schroder, GER 

Runge-Kutta Unvealed, J. Bohm, AUT 

The Horror Octahedron, W. Alvermann, GER 

RK6, Heinrich Ludwig, GER 

and others 

Impressum; 

Medieninhaber; DERIVE User Group, A-3042 Wiirmla, D'Lust 1, AUSTRIA 
Richtung: Fachzeitschrift 
Herausgeber: Mag. Josef Bohm 


The DERIVE-NEWSLETTER is the Bulle- 
tin of the DERIVE & CAS-TI User Group. 
It is published at least four times a year 
with a contents of 40 pages minimum. The 
goals of the DNL are to enable the ex- 
change of experiences made with DERIVE, 
TI-CAS and other CAS as well to create a 
group to discuss the possibilities of new 
methodical and didactical manners in 
teaching mathematics. 

Editor: Mag. Josef Bohm 
D'Lust 1, A-3042 Wiirmla 
Austria 

Phone: ++43-06604070480 

e-mail: nojo.boehm@pgv.at 
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Tania Koller, Vienna, Austria 
Dear Josef, 

Last week we wanted to produce a 3D-representation with DERIVE and four students 
found a very strange distorted view on their screens. It couid not be caused by VISTA 
because aii other screens (aiso VISTA) were ok. Fortunateiy I have some bright and ea- 
ger students, so I can send not only the riddle but also its solution which was found by 
one of my students in the weekend. Maybe that this could help other users, too. 
Regards 
Tania 

Correction of 3D-plots on Notebooks (Windows 7) 


eg. with a resolution of 1366 x 768 



1. Show Desktop - minimze all programs in order to show the desktop 
(shortcut WINDOWS + D) or rightclick on the right border of the systray - Show Desktop 



2. Rightclick on the Desktop - Resolution 







4. Click on „Keep Changes" (two black stripes appear on the screen). 



5. Open a new 3D-Plot Window - and now it should look fine! 


I|ii I 

E«t £« hMM Sm ]gM*> 

Oo^y#% s o0c>*;<^ Q--T 1 XX » 





JJiJd/JillJ iJiJ. 


Many thanks to Philipp Wietter, student of Tania Koller at HAK St. Polten. 
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Dog and Biker 

Roland Schroder, Celle, Germany 

There are dogs which become aggressive when they discover a biker. We will not investigate 
why this is the case. But we will discuss the question: what is the chance of a biker to escape 
an aggressive dog? We can imagine the situation as follows: 



The biker moves with constant velocity starting from point (-40,0) in direction to the origin and 
goes ahead in this direction to escape the dog if possible. The dog is watching his surround- 
ings in position (0,40) when he is discovering the biker at (-40,0). The dog is silly enough not 
to bar the biker's way but he undertakes 4m jumps (with also constant velocity) in the direc- 
tion of the - changing - position of the biker. If both have the same velocity the dog will never 
reach the biker because of his wrong strategy. 

If the dog is faster than the biker he will win but the dog gives up “hunting” the biker after 20 
jumps. What is the ratio of velocities giving the biker a chance to escape? 

We choose the time steps so that one jump of the dog takes one step. The velocity of the 
dog is the k-fold of the biker's velocity with /c > 1 to give the dog a chance. The positions of 
the dog and the biker at time t are illustrated by the following sketch: 

^ Dog [a,b] 


Biker [- 40 + 4(f - 1 )//c, 0]] 
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The position of the biker only changes in its x-coordinate, which is described by the DERIVE- 
expression 


#1: -40 + 4-t/k. 


The segment (the vector) from Dog to Biker is given by [#1,0] - [a,b] = [#1 - a, - b]. The 
length of this segment - this is one jump of the dog - should be 4 m. So we have to normal- 
ize the segment (= divide by its length) and multiply by 4. We use the factor: 

#2: 4/ABS([#1 - a, -b]) = #2: 4/ABS([-40 + 4-t/k - a, -b]). 

The jump starting from [a, b] is presented by #2 • [#1 - a, - b]. This jump starts at [a, b] and 
ends at [a, b] + #2 • [#1 - a, - b] = [a + #2 • (#1 - a), b - #2 b]. We need a third component in 
this vector, t + ^, which counts the jumps of the dog. We find by recursion the immediate 
successor element of [a, b, t] in [a + #2 • (#1 - a), b - #2 b, t + 1]. The ITERATES-command 
is the perfect tool to perform the recursion: 

#3: D(k):=ITERATES([a + #2 • (#1 - a), b -#2 b, t + 1], [a, b, t], [0, 40, 0], 20). 

For plotting we have to remove the time component. The sequence of points - landing points 
of the dog - are given by: 

#4: Dog(k) == VECTOR([(D(k))ini1, (D(k))ini2], n, 21). 

But we want to see the position of the biker, too. So we enter 
#5: Bike(k) == VECTOR([-40 + 4-(t- 1)/k, 0], t, 21). 

#6: Dog(1.2) 


#7: Bike(1.2) 

Now plot #6 and #7 (for k = 1 .2): 



D 


-N-L#75 


Roland Schroder: Dog and Biker 


P 7 


We can use a utility from DNL#62 to perform something like an animation - we can imagine 
following the jumps: 

moveptsCli'st) := VECTOR (list -IF(t = i , 1, ■») , DIM (list)) 

#8: i 

traceptsCli'st) := VECTOR (list -IFCt> 1, DIMCli'st)) 

#9: i 

#10: [movept5(Dog(l. 2)) , moveptsCBikeCl. 2)) ] 

#U_: [traceptsCDogCl. 2)) , traceptsCBi'keCl. 2)) ] 

Introduce a slider for t with 1<= t <= 21 (20 Intervals) and now you can either perform one 
jump after the other (#1 0) or follow the trace of dog and biker (#1 1 ). 


t = 10.00 


1.00 I 


1 



Dog 


30 

25 


20 


15 


10 


Biker 


45 -40 -35 -30 -25 -20 -15 -10 


5 


-5 5 


10 15 



I (Josef) wanted to transfer this nice problem to TI-NSpireCAS because recursion can be 
performed on a spreadsheet in an easy way. I wanted to have a very flexible model with vari- 
able positions of dog and biker and including a slider for the parameter k to investigate its 
influence on the success of the dog - or the biker. 


I came across one problem: the vector operations don't seem to work in the spreadsheet 
application (eg. norm(x,y) for finding the length of a vector). Philippe Fortin confirmed my 
conjecture and gave a hint to overcome this problem. I defined a function for the length of the 
vector and then the spreadsheet application accepted the function. 


You are invited to study the TI-Nspire-file. I started with creating the points for the dog and 
the biker and linking their coordinates to variables which I then used in the spreadsheet. 
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|| le [a, b,€,d)\ =nonn ( [c -a d-b ] ) 

Done 




m 

3xb 

Syb 

Bxd 

Fjyd 

♦ 







T 

1 ,200000,., 

1.000000,., 

-40.0000,., 

o.oooooo... 

o.oooooo... 

40.00000,., 

2 


2.000000,., 

-36.6666,., 

o.oooooo... 

-2,82842.,, 

37.17157,., 

3 


3.000000,., 

-33.3333,., 

o.oooooo... 

-5,52111.,, 

34.21363,., 

4 


4.000000,., 

-30.0000,., 

o.oooooo... 

-8,04423.,, 

31.10978,., 

5 


5.000000,., 

-26.6666,., 

o.oooooo... 

-10,3506.,, 

27.84170,., 

6 


6.000000,., 

-23.3333,., 

o.oooooo... 

-12,3730.,, 

24.39064,., 

7 


7.000000,., 

-20.0000,., 

o.oooooo... 

-14,0126.,, 

20.74208,., 

8 


8.000000,., 

-16.6666,., 

o.oooooo... 

-15,1219.,, 

16.89899,., 



DOG 


k := 1.6 
.5 2 . 


□ 

^ (12.8,34.6) 

□ 

□ 


(-30, o) 


□ 


5 


□ 


□ 


■o— *-9..— o- 


-•-o 

5 


-o— 

□ 
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Application of Computer Algebra (ACA) 2009 

Session: Applications and Libraries development in Derive 
25-28 June 2009, Ecole de technologie Superieure, Montreal, Quebec, Canada 


Another Look at a Trusted Mathematical Assistant 


Michel Beaudin 

Service des enseignements generaux, 
Ecole de technologie superieure (ETS), 
1100, Notre-Dame street west, 
Montreal, Quebec, Canada, H3C 1K3 
Email: michel.beaudin@etsmtl.ca 


Abstract 


• From the DERIVE user manual (version 3, September 1994), we ean read the 
following: “Making mathematics more exciting and enjoyable is the driving force 
behind the development of the DERIVE program”. In this talk, we will try to show 
how some mathematical concepts, studied by engineering students at university level - 
differential equations, multiple variable calculus, systems of non linear equations -, 
can be easily illustrated by DERIVE. Some will object that any other CAS could do 
the same: well, this is probably true but, according to us, not as quickly and naturally: 
“To accomplish this DERIVE not only has to be a tireless, powerful and 
knowledgeable mathematical assistant, it must be an easy, natural, and convenient 
tool”. Consequently, time can be spent to prove some theorem or formula and the 
computer algebra system helps to reinforce the mathematical concepts. Our examples 
will also make use of new features added in the latest version of DERIVE (version 
6.10 released in October 2004); features that were not exploited as should be - 
DERIVE has never been enough documented. But we are still convinced that Derive 6 
was “Far too good just for students” 

(http://www.scientific-computing.com/scwmarapr04derive6.html ). 
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1. Introduction 

This main goal of this paper is to show how DERIVE can be used at university level, for teach- 
ing to engineering students. For the past 15 years, we have attended many conferences about using 
CAS in teaching mathematics: when Derive was used in a presentation, it was — in the majority of 
cases — for high schools/colleges level presentation. This is normal because DERIVE, compared to 
“big” CAS like MAPLE or MATHEMATICA, is very simple to use and, at lower levels, this is impor- 
tant. But, at university level, even if DERIVE is a small system compared to the big ones, it can be 
used to explore many avenues, in domains such differential equations and multiple variable calculus — 
where, usually, one will use another system, more “powerful”. In fact, as far as undergraduate 
mathematics are concerned — the situation would probably be different at research levels — , DERIVE 
is so simple to use that if some mathematics teacher wants to go from the blackboard to the computer 
— without having to load an already prepared file — , DERIVE is an excellent choice: no waste of time 
by typing complicated or long commands in order to illustrate daily concepts. This is important for us: 
we are still teaching mathematics and we need software that will do the job correctly. But we also 
need software that can do heavy computations when this is needed. DERIVE has done the job since 
the last 15 years. We would like to be able to continue to use DERIVE for the upcoming years. Be- 
cause of its discontinuation, we are afraid that these good moments are behind us: only time will tell. 


2. The importance of symbolic integration of piecewise continuous functions 
in differential equations 


If we don’t pay attention to the endpoints of an interval [a, b\, then the indicator function, defined by 
CHI(a, X, b), will be the following function: 


Cm{a,x,b) 


[l if a <x<b 
[O elsewhere 


When you are teaching differential equations, you often need to deal with piecewise continuous func- 
tions, especially when you study Laplace transforms: each piece is defined on an interval, 2 distinct 
intervals don’t overlap, so CHI is very useful. So, where does the CHI function play an important 
role? To answer this question, let’s take a look at the following problems from a differential equations 
(DE) course. The first two ones are usually solved by taking the Laplace transform of both sides of 
the DE. The third one is done by computing the Fourier coefficient of the given signal or by using a 
table of Fourier series. The fourth one is not related with the CHI function but because it is usually 
solved by using the method of undetermined coefficients, this will become an opportunity to explore 
what DERIVE is doing when its “DSOLVE2” command is used. 
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a) Solve y" + 2y' + y = f (0? >"(0) - 0? - 0 where / is the function in figure 1. 



b) Thinking of a mass-spring problem, find the solution of 

y" + Ay = 505{t-7i)-\005{t-27i\ y(0) = l0, /(0) = 5, 

where S is the «Dirac delta function” (not defined in DERIVE). 

c) Find the Fourier series of the signal — square wave — shown in figure 2. What happens to the 
partial sums near a point of discontinuity? 


47T 


-Btt 


-2n 


2tt 


3tt 


At 


Figure 2 


d) Identify the steady-state solution in the mass-spring problem governed by the DE 


y* + ^y + ^9 y = 50co^{cot). 


Here co is some positive parameter. 
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Now, let us give the “DERIVE solutions” of each of these problems. For problem a), one first needs 
to note that the function in figure 1 is defined by 

0 if f < 0 

6t if 0 < f < 1 

<6 ifl<f<3 

24 - 6f if 3 < f < 4 
0 iff >4 


Using Laplace transform techniques, we can find the following solution: 

y{t) = g{t)u{t) - g{t - \)u{t - 1) - g{t - 2>)u{t - 3) + g{t - 4)w(f - 4), 

where g(f) = 12e"' +6fe”^ -12 + 6f, w(f) = STEP(f). Using DERIVE, the command 
“DSOLVE2 IV” does the job, because this command uses the method of variation of parameters in 
order to find a particular solution and the software has no problem to integrate piecewise continuous 
functions. So, we define / in line #1 (figure 3 below): this is the “input” and this is done with linear 
combination of CHI functions because the intervals don’t overlap. And use the “DSOLVE2 IV” com- 
mand to get the answer or “output” (no need to simplify this to get the graph if the option “simplify 
before plotting” is on). Of course, if one simplifies line #2, an answer involving SIGN functions will 
appear on the screen — because CHI functions are defined using STEP and STEP comes from SIGN!. 

#1: f ;z t, 1) + t, 3) + (24 - 6.t).}((3, t, 4) 

#2: DS0LVE2_IV(2 , 1, f, t, 0, 0, 0) 



Figure 3 
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For problem b), students learn that an expression of the form A5{t — t^^ means that an infinite force 

was applied at time t = t^ with total impulse of A (N.s). The expression AS{t-t^^ can be defined 

by the following limit, which does not exist in the classical way. 

A5it-ti^) = lim ^ — +^). 

^^0^ s 

But if one solves the differential equation with A5it — t^\ replaced by ^ — +6*) and 

s 

then let ^approach 0 from the right, he/she will get the right answer! Of course, the beauty of Laplace 
transform is that students don’t have to use that kind of trick. They simply use the fact that 
and find the solution. In order to fully appreciate the effect of the “Dirac delta func- 
tion”, the slider bar of DERIVE can be used if one wants to understand why the solution (line #5 in 
figure 4) is obtained by taking some limits (this has to be done “live” during a talk). The graph of the 
solution (line #5) is also shown after line #5 (note: there is no discontinuity!). 

50 100 

#1: input := "ti tt + a) - t, Z-jt + b) 

a b 

#Z: output := DSOLVEZ_IV(0 , 4, input, t, 0, 10, 5) 

rea1_ouput ;= 1im 1im output 
b-s-0+ a-s-0+ 

#4: rea1_ouput 

Z5.SIGN(t - 7T).SINCZ.t) 

#5: - Z5-SIGNCt - Z -tt) ■ SIN(Z ■ t) + + lO-COS(Z.t) - 10-SINCZ-t) 

Z 



Figure 4 
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For problem c), because DERIVE has a “Fourier series” built-in function, it is quite easy to obtain 
partial Fourier sums. First, we define the function on its interval of definition, using CHI again and 
then, extend it periodically using the modulo MOD function of DERIVE. Let us recall (fact learned at 
the first international DERIVE conference, in Plymouth, UK, in 1994!!!) that if g(x) is a defined func- 
tion of the variable x, defined over the interval [a, b\ — using possibly one or many indicator functions 
— , then we can extend g periodically over the entire real line, period being b - a, with the help of the 
definitions^) = g(mod(x - a, b - a) + a). A simple partial sum (line #3 and its simplification, line #4) 
tells us that the partial sum of order 2n - I is given by an expression of the form #5. Differential cal- 
culus tells us that the first critical point after 0 is located at x = n/li^ri), it is a local maximum — second 
derivative test shows this — whose y-coordinate is about 5.35795949. DERIVE won’t find the limit, as 
n goes to infinity, of expression #5 when x = nil {In). This can be done “by hand”, on the blackboard, 
and we can show that this limit is exactly 3 + 4/;rSI(7r), where SI stands for Sine Integral. 



For problem d), we said that DERIVE is using, within its “DSOLVE2” command, the method of varia- 
tion of parameters. This gives a very ugly answer compared to the one we get if we use the method of 

undetermined coefficients. In fact, a particular solution of y" + + ^9 y = 50co^{^cot) will be of 

the form a cos(gj t) + p sin(G7 1). Simple computations — but we did it using DERIVEl — show that 
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a) How can we get DERIVE to obtain one complex solution of an equation like the following 
one: 1" +21 Those who are familiar with the LambertW function of Maple know that 

there are many complex solutions for an equation of the type a"" . 


b) According to the theorem about existence and uniqueness of a solution of a first order ODE, 

dy x^ + e ^ 

there is a unique solution to the problem — = , y(0) = 2, defined in some interval 

dx 6 + 2y 

about 0. How is this solution related with the implicit curve obtained by solving the ODE as a 
separable one? And how does the Euler numerical method compare? 

c) In DERIVE, level curves of a function of 2 independent variables have to be plotted in 2D, not 
on the surface in 3D. How can we visualize this in 3D, using the slider bar? 



Figure 6 


But if we are searching a complex solution, one needs to substitute for x the expression x + iy, and 
plot, in the same window, the implicit curves defined by RE(eq(x + iy)) = 0 and IM(eq(x + iy)) = 0 
where eq(z) = 2^ -2. The fast 2D plotter and Newton’s method (for 2 variables) will give any 

complex solution (or Newton’s method for 1 variable, starting point being complex). Finding a com- 
plex solution using the real and imaginary parts of the given equation is a nice way to use complex 
numbers (" w e C = 0 RE(w) = 0 a IM(w) = 0"). It would have been a good addition to “version 7 of 
DERIVE'' to have, within the solve command, the possibility, for a system of 2 equations, to use a 
starting point without calling the “Newtons” function. 
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For problem b), let us note first that we can find a neighbourhood of the point (0, 2) — obviously lim- 
ited by the horizontal line y = -3 where the derivative becomes infinite — where both the function 

and the partial derivative with respect to y are continuous. So, there exists an interval about 

6 + 2y 

the point 0 and a unique function that satisfies the ODE and the initial condition. The “separable” 
function of DERIVE is a good choice to obtain the solution, defined implicitly by 

x^ 

y^ + 6y = -e + 17. Now, the importance of a 2D plot window, where all kind of plots are possi- 

ble, becomes clear. Solving for y — of course, here it is possible with the well-known second degree 
formula — and keeping the “good” branch leads to 

- x/2 3 x/2 

' CyCe • - 3) - 3-^3 -e ) 

cpCx) := 

3 

Now, in the same window, here are plots of the implicit solution, the explicit solution (function 
^(x)and points generated by the Euler’s numerical method: it is always good to be able to see the 
numerical solution (Euler’s method in this case) and the exact solution (implicit and explicit when 
possible) in the same window. 



Figure 7 
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For problem c), let us recall that if f = fix, y) is a scalar field, then the (implicit) curve defined by 
fix, y) = c is called a level curve of / (c being the level). The “old” feature of DERIVE that consists to 
split the screens can be used to see a particular level curve in 2D and in 3D: for the 3D view, simply 
intersect the surface z = fix, y) with the plane z = c. Here is an example with the function 

f(x,y) = x~— Let’s plot the level c = I curve, in 2D and visualize it in 3D (note that one has 

to make the box turn if axes are to be oriented like they are in 2D): 



Figure 8 


By the way, beginning with DERIVE 6, a Grobner basis function (groebner basis) started to be used 
by the author when heavy polynomial systems needed to be solved. That became a nice opportunity to 
extend the famous row-reduce function used with linear systems. 


4. 3D plotting facilities: quite acceptable 


Before version 5 of DERIVE, 3D plots were very limited. For instance, figure 8 shows, on the 
right, something that it was impossible to do. Hopefully, David Parker was hired by Texas Instru- 
ments to allow DERIVE 3D plotting capabilities to be compared with ones of big major CAS. If 
DERIVE 6 does not support implicit 3D plotting, it has parametric 3D plotting — one parameter for a 
space curve, two parameters for a surface — that does a good and correct job. Let us concentrate over 
a single, simple example, an example where DERIVE helps the mathematics teacher to concentrate on 
concepts and not on complicated commands. Suppose you want to generate a torus, taking the circle 
of radius 1, centered at the point (2, 0, 0) in the plane y = 0, and letting it turn around the z-axis. You 
should get something like this: 
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Figure 9 


In DERIVE, plotting the expression (vector) [2 + cos t, 0, sin t\, where t ranges from 0 to 
Ik, will give the circle — of course, one has to find parametric equations for the circle, but students 
learn how to do this in 2D, so it is quite easy, here! Using the matrix “ROTATE Z(^)” will finish the 
job. Note that we had to transpose the matrix obtained by the vector in #1 of figure 10 in order to 
perform the matrices product in #2 and, then, we had to put the result in vector form in order to get the 
3D surface. Before doing this, one can use the slider bar in #5 of figure 10— replacing shy a for ex- 
ample — and see the circle turning around the z-axis: this is something that has to be done “live” with 
DERIVE and if you want to do it, you will see that you don’t have to type too many commands. With 
DERIVE, it is easy and fast to get this. Figure 10 concludes this example. 
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#1: [2 + COS(t), 0, SIN(t)] 

z 

A 

X y 



# 2 : 


#3: 


#4: 


# 5 : 


R0TATE^(s).[[2 + COS(t), 0, SIN(t)]]’ 

■ C0S(s).(C0S(t) + 2) ■ 

SIN(s).(COSCt) + 2) 

SIN(t) 

■ C0S(s).(C0SCt) + 2) ■ 

SIN(s).(COSCt) + 2) ’ 

SIN(t) 

[[COSCs).(COSCt) + 2), SIN(s).CCOSCt) + 2), SIN(t)]] 

: z ; 

: A I 


X y z 



Figure 10 
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5. Conclusion 

This additional look at a trusted mathematical assistant — limited look where examples con- 
cern mainly mathematics at university level — is my way to say, again, thank you to Albert Rich and 
David Stoutemyer, the “fathers” of DERIVE. Also, I want to thank Theresa Shelby, for her nice job 
regarding the windows versions that appeared after version 3. These persons have created a piece of 
software that gave me so much enthusiasm during the past 18 years of teaching mathematics to engi- 
neering students at ETS, in Montreal, Canada. Also, the TI-92 calculator — now Voyage 200 — is, as 
far as we are concerned, a direct result of DERIVE' s spirit (I know that Albert Rich won’t agree!). 
And, since 1996, we have started to use both (see [3]): that became more interesting in 1999 when 
every student at ETS started to buy — compulsory purchase — a TI-92 Plus (Voyage 200 from 2002) 
because, as a math teacher, we were able to switch from the calculator to DERIVE easily. Also, for 
many problems, the calculator was sufficient: this was a good reason to continue to think that a laptop 
computer is not necessary for a student when attending a mathematics course. 

DERIVE is no longer on the market. A certain part of its spirit is somewhere in Voyage 200 
calculator, whose future is uncertain... Nspire CAS software pretends to become the true successor of 
DERIVE. Again, only time will tell. 
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Abstract 

This paper establishes the theoretical aspects which have been considered in order to elabo- 
rate the Random_distiributions package for Derive 6 as well as the description of the different 
algorithm developed in the package. In section the theory on random number generation is 
presented (from |Rnbinstein, 1981] ). After explaining Derive’s random function (section o> 
the more efficient algorithms ran2 and mzranlS are developed (section [L^ and [L3| respectively) . 
Section presents three different general methods for generating continuous distributions to- 
gether with one for generating discrete distributions. Section is dedicated to describe different 
algorithms for generating random values from continuous distributions (Uniform, Exponential, 
Normal, Lognormal, Weibul, Gamma, Beta, Chi-square, Student’s t, F, Z, Pareto, Logistic, 
Cauchy and Irwin-Hah distributions). Section presents different algorithms for generating 
discrete distributions (Uniform discrete, Bernouille, Rademacher, Binomial, Poisson, Geomet- 
ric, Negative Binomial and Hypergeometric distributions). In section some algorithms for 
generating different distributions by approximations are developed. Section is devoted to de- 
velop some graphical approaches in order to check graphically how the generated sample fix the 
distribution. Finally, in sections and some examples and possible extensions to this work 
are shown. 


1 Random Number Generation 

The most commonly used methods for generating pseudorandom numbers are congruential genera- 
tors. A congruential method is one that produces a nonrandom sequence of numbers according to 
some recursive formula based on calculating the residues modulo some integer m of a linear transfor- 
mation. It is readily seen from this definition that each term of the sequence is available in advance, 
before the sequence is actually generated. Although these processes are completely deterministic, it 
can be shown that the numbers generated by the sequence appear to be uniformly distributed and 
statistically independent. 

Congruential methods are based on a fundamental congruence relationship, which may be ex- 
pressed as: 

Xi^i = aXi~\-c (mod m) i = 0,1,2, ... ,n (1) 

where the multiplier a, the increment c and the modulus m are nonnegative integers. 
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Given an initial starting valne Xq (called the seed), 0 yields a congrnence relationship (modulo 
m ) for any value i of the sequence {Xi}. Generators that produce random numbers according 
to Q are called mixed congruential generators. The random numbers on the unit interval (0, 1) 
can be obtained by: 

Ui = ^ ( 2 ) 

m 

Glearly, such a sequence will repeat itself in at most m steps, and will therefore be periodic. 

As Xi < m for all i, the period of the generator cannot exceed m, that is, the sequence A* 
contains at most m different numbers. Because of the deterministic character of the sequence, the 
entire sequence recurs as soon as any number is repeated. The sequence is said to get into a loop, 
that is, there is a cycle of numbers that is repeated endlessly. Modulus m should be chosen as large 
as possible and appropriated values of a and c in order to make the period p maximum (that is, 
p = m) must be found. When this happens the random number generator has a full period. It can 
be shown that the generator defined in ([^ has a full period m, if and only if: 

1. c is relative prime to m, that is, c and m have no common divisor. 

2. a = 1 (mod g) for every prime factor g of m. 

3. a = 1 (mod 4) if m is a multiple of 4. 

Since most computers utilize either a binary or a decimal digit system, the best selection for m 
is m = 2^ or m = 10^, respectively where f3 is the word-length of the particular binary or 
decimal computer. 

For binary computers, in order to develop a full period generator when m = 2^, the parameter 
c must be odd and a = 4fc + 1 for some k gN. 

The second widely used generator is the multiplicative generator: 

Aj+i = aAj (mod m) z = 0, 1, 2, . . . , n (3) 

which is a particular case of the mixed generator 0 with c = 0. 

Another common type of generator in which depends on more than one of the preceding 

value^ For example: 

Aj_|_i = Ui Xi-j-^ + 02 Xi-j^ + • • • + Ufc ^i-jk T ^ (mod m) or 

Ai+i = a Xi_j^ ■ Xi_j^ ■ ■ ■ Xi_j^^ + c (mod m) 


Nowadays “the best” generators use combinations of the generators described along this section 
in order to increase the randomness and the period of the generated sequences. 

^ These generators are often called Fibonacci generators because one example is given by the Fibonacci serie: 

Xi+i =Xi+ (mod m) 
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1.1 Derive’s random function 

Derive’s random function uses a mixed generator given by: 

Xi+i = 2, 654, 435, 721 + 1 (mod 2^^) 

which satisfies the conditions to be a full period generator, that is, the period of Derive’s random 
function is 2^^ = 4, 294, 967, 296. 

Derive’s random function RANDOM (n) can be used with any n G Z with the following meanings: 

• If n > 1, RANDOM (n) simplifies to a random integer in the interval [0,n). 

• RANDOM(l) simplifies to a random number in the interval [0, 1). 

• If n < 0, RANDOM (n) simplifies to —n and initializes the random number state variable to 
—n. 

• RANDOM (0) simplifies to the time in centiseconds since the current calendar year began and 
initializes the random number state variable to that time. 

Although this is a “good” generator, the following two subsections describe two different algo- 
rithms, implemented in the package Random_distributions, which periods are quite much longer 
and also improve the randomness. 


1.2 ran2 algorithm 


The ran2 algorithm was proposed by L’Ecuyer and is described in [Press and Teukolsky, 199~2] 
Press et ah, 1999| . 


and 


This algorithm merges the following two multiplicative generators: 


Xj+i = 40014 Xi (mod 2^^ -85) 
Yi+i = 40692 Yi (mod 2^^ - 249) 


This algorithm has been used for a long time as one of the best generators and its period is 
about 2.3 • 10^® = 2, 300, 000, 000, 000, 000, 000 which is more than 535, 510, 480 times longer than 
Derive’s random generator period. 

The implementation on Derive has been carried out “translating” the “C” code developed in 
Press et ah, 1999| and it uses the following two functions: 

• ran2(n) which is the main algorithm. This function returns a vector of size n of random 
numbers in the interval [0, 1). 

• ran2_initialize() This auxiliar function is used to set the variables and constants needed 
for the algorithm. 


1.3 mzranlS algorithm 

The mzranlS algorithm was proposed by G. Marsaglia and A. Zaman as an alternative to ran2. This 
algorithm is described in [Marsaglia and Zaman, 1994| . 

This algorithm merges the two generators: a mixed one with a Fibonacci’s like one. 

Ai+i = 69069 Ai + 1,013,904,243 (mod 2^^) 

F,+i = - Ti _2 - “c” (mod 2^2 _ ig) 

where the second one is a subtract- with-borrow generator (because of the term “c”). 
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This algorithm has been found to be at least as good as ran2 but simpler, much faster and with 
periods “millions and millions” of times longer. Specifically, its period is over 

2^4 = 19^ 807, 040, 628, 566, 084, 398, 385, 987, 584 


that is, 8.611.756.795 times longer than ran2’s period and 4,611,686,018,427,387,904 times 
longer than Derive’s period. 


The implementation on Derive has been carried out “translating” the “C” code developed in 
IMarsaglia and Zaman, 19M| and it uses the following function: 


• mzranl3(n). This function returns a vector of size n of random numbers in the interval 
[0,1). Previously, when the package Random_distribution is loaded, the needed constants 
and variables are initialized. 


This algorithm is the base for all the random distribution generations developed in this package. 
This, in the following, when we say that a value u is generated from 7/(0, 1), this is done using 
mzranlS algorithm. 


2 Different methods for random variate generation 

This section presents some general methods for generating random variables from different contin- 
uous and discrete distributions. In the following subsections three general methods for continuous 
distributions and one for discrete distributions are described. 

2.1 Inverse transform method 

Let T be a random variable with cumulative probability distribution function Fx {x) . Since Fx {x) 
is a nondecreasing function, the inverse function Fx^ {y) may be defined for any value of y between 
0 and 1 as: Fx^ {y) is the smallest x satisfying Fx{x) > y, that is, 

Fx^{y) = inf {Fx{x) >y} , 0<y<l 

If U is uniformly distributed over the interval (0, 1), then X ilA). So, to get a value x 

of the random variable T, a value u from a random uniform variable 7/(0, 1) can be obtained 
and compute x = Fx^ (y). Thus, the general algorithm for the inverse transform method is: 

1. Generate a value u from 7/(0, 1). 

2. Obtain x = F^^{u) as the random number from the variable X. 

The only condition needed for this method is that Fx' exists in an analytical form. 

The following Derive’s program has been developed in the package Random_distributions to 
obtain a formula to generate X using the inverse transform method: 


Inverse.transf orm_method(f , ini := 0, u) : = 
Prog 
( 

u : G Real [0,1], 

Solve (u = INT(f, X, ini, x) , x. Real) 

) 
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2.2 Composition method 

This method is employed by Butler and consists of expressing the probability density function fx (x) 
of the distribution to be simulated as a probability mixture of properly selected density functions. 

Let g(xji/) be a family of one-parameter density functions, where y is the parameter identifying 
a unique g{x). If a value of y is drawn from a continuous cumulative function Fy{y) and then if 
X is sampled from the g{x) for that chosen y, the density function for X will be 

fx{x) = j g{x\y) dFy{y) 

If y is an integer parameter, then 

fx{x) = ^Pig{x\y = i) 

i 

where 

= 1 ; P,>0 ; P, = P[y = i] i = l,2,... 

i 

This method may be applied for generating complex distributions from simpler distributions 
that are themselves easily generated by the inverse transform method or by the acceptance-rejection 
method described below. 

2.3 Acceptance— rejection method 

This method is due to von Neumann and consists on sampling a random variate from an appropriate 
distribution and subjecting it to a test to determine whether or not it will be acceptable for use. 

To carry out this method, the probability density function fx{x) of the variable X must be 
expressed as: 

fx{x) = C -h{x) -g{x) 

where C >1^ h{x) is also a probability density function, and 0 < g{x) < 1. After generating two 
random values u and y from W(0, 1) and h{y), respectively, the test to see wether or not the 
inequality u < g(y) holds must be done, and: 

1. If the inequality holds, then accept y as a, variate generated from fx{x). 

2. If the inequality is violated, reject the pair u , y and try again. 

So, the general algorithm for the acceptance-rejection method is: 

1. Generate a value u from W(0, 1) 

2. Generate y from the probability density function h{y). 

3. If u < g{y), return y as the variate generated from fx{x) 

4. Go to step 1. 
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2.4 Inverse transform method for discrete distributions 

The inverse transform method is the easier method to use not only for continuous distributions but 
also for discrete distribution. 

Let T be a random discrete variate which finite or infinite possible values are 

, X2 1 • • • 1 •)•••• 

Let Fx{x) be its probability mass function given by 

F{xi) = P[X = Xi] = Pi i = 

The inverse transform method can be described as follow: 

1. Generate u from 

2 . i := 1 . 

3. p :=pi. 

4. If u < p deliver Xi as the generated value. 

5 . i:=i + l. 

6. p := p + pi 

7. Go to step 4. 


On the other hand, the values Xi can be assumed to be all integers and Xj+i = a;* + 1, because 
if this is not the case, the correspondence 4>{xi) = i can be established and consider a new random 
discrete variate y which values are 1,2,... and Fy(i) = P[y = i] = P[X = Xj] = Pi, which is 
equivalent to X and verify the above condition. 

Let X = {ini, ini + l,ini + 2,...} for some ini G Z with probability mass function 
Fx{x) := P[X = x\= Px ; x = ini, ini + 1, ini + 2, . . . 

The following Derive’s program has been developed in the package Random_distributions to 
generate an element of X using the inverse transform method: 

random_discrete_aux(F, ini := 0, aleat, i_, p) : = 

Prog 

( 

i_ := ini, 

p := SUBSTCF, x, ij , 

Loop 

( 

If aleat > p, 

RETURN i_), 
i_ := i_ + 1, 
p := p + SUBSTCF, x, iJ 

) 

) 


while the following Derive’s program generate a sample of size n of X. 
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random_discrete(n ;= 1, F, ini ;= 0, vecaleat) ;= 

Prog 

( 

vecaleat := random_uniform(n) , 

VECTOR(random_discrete_aux(F, ini, vecaleat sub j), j, n) 

) 


3 Continuous distributions random generation 

This section describes generating procedures for different continuous distributions. 

All these continuous distributions are presented with their probability density functions and 
with a drawing of a bars diagram, obtained by the Derive’s algorithm developed in the package 
Random_distributions, together with the plot of the corresponding probability density function in 
order to see graphically if it is a “good” sample of generated values (see section 

See IRubinstein, 198l| , |Galan, 1991| and [Wikipedia, 2009| for further information on the algo- 
rithm described below. 

3.1 Uniform distribution 


A random variable X has an uniform distribution in 
the interval (a, 6) ( X U (a, b) ) if its probability 

density function is: 


1 


fx{x) = 


b — a 

0 


X e (o, h ) 

otherwise 



W(0,1) 


To generate X, first u must be generated from ZY(0, 1) and then return 
the generated value. 


a + {b — a)u 


as 


3.2 Exponential distribution 



A random variable X has an exponential distribution 
with parameter A > 0 ( A T(A) ) if its probability 
density function is: 


fx{x) 


\ e“^/^ a: e [0, oo) 
A 

0 otherwise 


To generate X the inverse transform method is used: 

U = Fx{x)= f \ e“*/^ dt = 1 - e“^/^ ^ X = -X \n{l - U) = -X In (U) 

Jo A 


Thus, to generate X, u must be generated from W(0, 1) and then return 
the generated value. 


X = —X Inn 


as 
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3.3 Normal distribution 


A random variable X has a normal distribution with 
parameters /r and cr ( A ^ A/"(/i, cr) ) if its probability 
density function is: 


fx{x) = 


1 




e 2 <t 2 


X e 



A^(0, 1) 


To generate A, Box and Muller theorem establishes that: 

If and U2^U{Q,1) then 

Zi = y/— 21n(^Yi) cos(27tZY 2) and Z 2 = y/— 2 In(^Yi) sin(27rZY2) 
are independent standard normal deviates A/”(0, 1) 


On the other hand, ii Z A/"(0, 1) then a Z) 

Thus, to generate A, u\ and U 2 must be generated from U{0, 1) and then return any of 
the values /i + ay/— 21n(ui) cos( 27 t U2) or /i + ay/— 21n(ui) sin(27r U 2 ) as the generated 


value. 

3.4 Lognormal distribution 



0.p^9 

/ M 


C|'.374.4 


1 


b.li>498 


0. 



8 1.6 2.4 3.2 


4.8 5.6 6.4 7.2 8 


£AA(0, 1) 


If the random variable Z cr) then A = e^ 

has the lognormal distribution with parameters fi and 
a ( A ). Its probability density function 

is: 


1 


fx{x) 


xa y/^ 

0 


(In x — iJb)^ 

e 2^2 X e [0, oo) 


Otherwise 


To generate A, 2 : must be generated from A/”(/i, a) and then return 
generated value. 


as the 


29 


3.5 Weibul distribution 


A random variable X has a Weibul distribution with 
parameters a > 0 and > 0 { X ^ W{a,(3) ) if its 
probability density function is: 


fx{x) 


^ ^ Q X G [0, oo) 

jja L > / 

0 Otherwise 



To generate X the inverse transform method can be used: 


U = F^{x)= ^ e-(^ dt = 1 - e-(f)“ 

^ ^ Jo /?“ 


v\ « 

-J = -ln(l-W) = -lnW'-^T(l) 


X = (3 


1/a 


To generate X, v must be generated from T(l) and then return x = 
generated value. 


as the 


3.6 Gamma distribution 



A random variable X has a Gamma distribution with 
parameters a > 0 and P > 0 { X ^ P) ) ii its 
probability density function is: 

{ r^a 1 ^ x/ (3 

— 

0 


ay J 


Otherwise 


The inverse transform method cannot be applied since F;^^{x) does not exist in an explicit form. 

In this case, different algorithms have been developed in order to generate samples of the Gamma 
distribution. Finally, the main algorithm chooses which is the appropriated one depending on the 
values of parameters a and p. 

3.6.1 random -gammal 

This algorithm is valid for values of a > 1. The following two properties of gamma distribution are 
the base to develop this algorithm: 

1 . g(l,l3)=S{l3). 

2. If Xi Q{ai,P) and X 2 Q{a 2 ,P) then X = Xi + X 2 ^(oi + 0 : 2 ,/?), that is, gamma 
distribution is reproductive with respect its first parameter. 
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Let a > 1 , m = floor(o;) and S = a — m where floor(a) is the integer part of a. In order to 
generate X^Q{a,(3), a mixture of Q{m,f3) and Q{m + 1,(3) with probabilities 1 — 5 and 5 
respectively can be used. On the other hand, in order to generate them, m or m + 1 variables from 


Q{1^(3) = S{I3) must be generated as shown in section 3.2 Thus, given ■ ■ ■ itim+i W(0,1): 

T = -/31n(^Yi) -/31n(^Y2) - 


/31n(Z4^) = -/31n j ^ ^(m,/3) and 

( m+1 


Q{m + l,/3) 


The algorithm is: 

1. Get «, til, • • • , 'Wm, 'Wm+i from ^Y(0, 1). 


2. Let X = Ui- 


i=l 


3. If 5 < M then let x = x • Um+i- 

4. Return — /5 ln(a:) as the generate value. 

3.6.2 random _gamma2 

If 0 < a < 1 then X = y-V where X^Q{a,f3) ; y'^Be(a,l — a) and V^S(/3). 

(Beta distribution Be is described in section |3?7| ) . Thus, random_gainma2 algorithm to generate 
T Q{a,(3) (0 < a < 1) can be described by: 


1. Generate y from Be{a, 1 — a) (using algorithm random_beta4 described in section 3.7.2). 

2. Generate v from S{l3). 

3. Return y ■ v as generated value. 

3.6.3 random_gamma5 

This is an acceptance-rejection method due to Gheng and describes gamma generation Q{a, 1) for 
a > 1. Let us remember that the acceptance-rejection method is based in the following decomposi- 
tion: 

fx{x) = C ■h{x) ■g{x) 

Gheng’s procedure uses: 


h{x) = 


C = 


g{x) = 


X = 


( ^ 

fl x^ ^ 

J (h 


1 0 


4a“ 

r(a) 

e“ A 

^a-X 

(/r -|- X' 

x/2a 

- 1 


x > 0 
Otherwise 


.2 e"' 


4 a“+^ 


where 


fi = a 


X 


Setting a = — ,6 = a — ln4 and c = a -|- a Gheng’s algorithm can be written as: 
A 
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1. Get u\ and «2 from fY(0, !)• 


Ml 


2. Let M = a In , , . 

\ 1 — Ml / 

3. Let X = a . 

4. If 6 + c M — a: > ln(MfM 2 ) retnrn x as the generated valne for Q{a,l). 

5. Go to step 1. 


3.6.4 random _gamma9 

This algorithm uses the approximation to Gamma distribution by Normal distribution when a and 
(3 are not “small” . 


a 


Specifically, if Z ^ M \\I^ \ — ] , \ — 1 then Q[a,(3) 


P 


1 1 


2a V a 


Thus, the algorithm 


is: 


1. Generate from A/" ( In ( — ^ — — , y — 


2. Deliver x = as the generated value for Q{a,P). 


3.6.5 random_gamma 

Finally, the following algorithm runs the above algorithms depending on the parameters a and P 
in order to generate a sample of size n from Q{a,P): 


random_gainma(n :=1, a := 1, P :=1) : = 

Prog 

( 

If (a = 1, return random_exponential(n, P)) , 

If (a < 1, return random_ganfflia2 (n , a, P)) , 

If {p = 1, return random_gamma5(n, a, /?)), 

If (a > 20 and (3 > 20 , return random_gamma9(n, a, P)) , 
return random_gainmal (n , a, P) 

) 


3.7 Beta distribution 


X has a Beta distribution with parameters a > 0 
and P > 0 ( T i3e(a, /3) ) if its probability density 
function is: 


fx{x) 


r{a + p) 
T(a) • T{i3) 


x^-\l-xf-^ 


X G [0, 1] 


0 


Otherwise 



Be{2,5) 


The inverse transform method cannot be applied since F^^{x) does not exist in an explicit form. 
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In this case, as in the case of Gamma distribution, different algorithms have been developed in 
order to generate samples of the Beta distribution. Finally, the main algorithm chooses which is the 
appropriated one depending on the values of parameters a and jS. 


3.7.1 random_betal 

This algorithm is based in the following result: 

If yi ^ Q{a, 1) and 3^2 1) then X = Be{a, (5). 

Ti + T2 

Its implementation is therefore trivial: 

1. Generate yi and r /2 from ^(a, 1) and respectively. 

2. Deliver x = — ^ — as the generated value for Beia,(3). 

Vl +V2 


3.7.2 random_beta4 


This algorithm has been developed for using in algorithm random_ganuna2. It is due to Johnk and is 
based on the following result: 


Let Ui W(0, 1) and U 2 7/(0, 1) and let = 7/^“ and 3^2 = T/a^^- If 3^i + 3^2 < 1 

then X = Be(a, (5). 

3^1 + 3^2 ^ ’ 

The algorithm is: 


U/2 liuiii 


1 . Generate u\ and 

2. Set U\ = and y 2 = U 2 

3. If yi + ?/2 < 1 deliver x = 

4. Go to step 1. 


1//3 


yi 


yi + 


as the generated value for Be{a,P). 


3.7.3 random_beta7 

This algorithm uses the approximation to Beta distribution by Normal distribution when a and (3 
are not “large enough” . 

X \ . ^ ^ a\ a — (3 


Specifically, if X then In 


1 - 


A/"(/i, (j) where = In — + 


/? J 2af3 


and 


G = 


la-\- (3 
af3 


Thus, 


In 


X 


1 — X 


= z 


X = 


(1 + e^) 


and hence, the algorithm is: 


1. Generate z from A/" ( In ( ^ ^ -i / 

a[3 


(3 J 2a(3 


2. Deliver x = 


(1 + e^) 


as the generated value for Be(a,(3). 


33 




3.7.4 random_beta 


Finally, the following algorithm runs the above algorithms depending on the parameters a and (3 
in order to generate Be(a,(3): 


random_beta(n := 1 , a 

II 

II 

II 

If (o;<4or/3<4, 


random_betal (n, a. 

(3), 

random_beta7 (n , a, 

P) 

) 



3.8 Chi-Square distribution 



Let 2i ^ 1) z = l,...,/c k standard normal in- 

k 

dependent distributions. In this case, A" = E has 


i=l 


the chi-square distribution with k degrees of freedom 

(X^x\k)). 


X^(35) 


Its probability density function is given by: 

^/ c / 2-1 q -^12 

fx{x) 


2k/2 Y (I) 


X e [0, oo) 


0 


Otherwise 


Although the algorithm to generate A X^ik) would be trivial by definition, it would need k 
values from A/"(0, 1) which require many operations if k is “large”. Thus, in the next two sections, 
two different algorithms which improve (in number of operations) the “definition algorithm” are 
described. 

3.8.1 random_chi_square2 

This algorithm is valid for “large” values of k (say k > 30) and it uses the following approximation 
from the standard normal distribution: 

If X x^{k) then Z = \J2X — \/2k — 1 is such that Z^M{0^1) 

{Z + ^/2k-l) 


Solving for A, A = 


from A/"(0, 1) and then return 


3.8.2 random_chi_square3 



Thus, to generate A X^{k) z must be generated 
as the generated value. 


This algorithm is based in the fact that x^(fc) is a particular case of a gamma density. Specifically, 
y^{k) = Q 2^ . Thus, to generate A 
return x = g I as the generated value. 


^^{k) = Q V Thus, to generate A X^{k), g must be generated from Q then 
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3.8.3 random_chi_square 


Finally, the following algorithm runs the above algorithms depending on the parameter k in order 
to generate 


random_chi_square(n ;= 1, k ;= 1) ; = 
If ( k > 30, 

random_chi_square2(n, k) , 
random_chi_square3(n, k) 

) 


3.9 Student’s t distribution 


Let Z A/"(0, 1) and y X^ik) independents. 


Then T = 




VyJk 

k degrees of freedom { X ^ t{k)). 


has a Student’s t distribution with 


Its probability density function is: 



-3.8 -3.04 -2.28 -1.52 -0.76 


0.76 1.52 2.28 3.04 3.8 


fx{x) = 


r(¥) 


1 + 


X 


2\ -(fc+l)/2 


X G 


t(35) 


v^r(l) ' k 

To generate X ^ t{k), z must be generated from A/"(0, 1) and y from y^{k) and then 


return 


X = 




as the generated value. 


3.10 F distribution 

Let Ti X^(^i) and T 2 X^(^ 2 ) independents. 

Then X = ^ has a T distribution with k\ and 

^ 2/^2 

k 2 degrees of freedom { X ^ T{ki,k 2 ) )■ 

Its probability density function is given by: 

^ ^ .. xG(0,oo) 

^(32,45) nW= r(f)r(f) + 

0 Otherwise 

To generate X ^ F(ki,k 2 ), yi and y 2 must be generated from and X^{k 2 ) 



0.14 0.28 0.42 0.56 0.7 0.84 0.98 1.12 1.26 1.4 1.54 1.68 1.82 1.96 2.1 2.242.38 2.52 2.66 2.8 2.94 3.0 


respectively and then return 


X = 


yi/ki 

1/2 7^2 


as the generated value. 
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3.11 Z distribution 


Let y y~{ki,k 2 ). Then, the variable X = 


In (3^) 


has a Z distribution with k\ and /c 2 degrees of 
freedom ( X Z(ki,k 2 ) ). 

Its probability density function is given by: 

kil2 


fx{x) = 


I) efa- 


X G 



X = 


To generate X ^ Z{ki,k 2 ), y must be generated from y{ki,k 2 ) and then return 
as the generated value. 


ln(y) 


3.12 Pareto distribution 



Pa(8, 1) 


A random variable X has a Pareto distribution with 
parameters a > 0 and a:o > 0 { X ^ Pa(«, Xq) ) if 
its probability density function is: 


fx(x) 


a /xq\“+^ 
Xq \ X J 


X e [Xo, oo) 


0 


otherwise 


To generate X, the inverse transform method can be used since: 

^ A = 




'xo 


Xo 


Xo 


u 


Thus, to generate X ^ Pa(a,a:o), u must be generated from U{0,1) and then deliver 
as the generated value. 


Xo 


l/a 


3.13 Logistic distribution 


A random variable X has a Logistic distribution with 
parameter a > 0 { X ^ kl{a) ) if its probability 

density function is: 


fx(x) 


a e 

(1 + 


X G M 



£( 10 ) 


36 


To generate T, the inverse transform method can be nsed since: 


U = FAx) = 


(x) = r 

j — 


a e 


-t 


dt = 


T = - In 


1 - 


/-inf (l + e-*)“+^ (l + e-^)“ V 

Thns, in order to generate X ^ T(o;), u must be generated from fY(0, 1) and then deliver 
as the generated value. 



3.14 Cauchy distribution 



A random variable X has a Cauchy distribution with 
parameters a > 0 and /5 > 0 { X ^ C{a,(3) ) if hs 
probability density function is: 


fx{x) = 




7T + (x — Ct)^] 


X G 


To generate the inverse transform method can be nsed since: 


U = FAx) = 




Loo vr [/?2 + (t - ay] 


atan ( 

21 dt = - + 


7T 


X = aX (3 tan 


TT [U — 


Thus, in order to generate X ^ C{a, P), xl must be generated from W(0, 1) and then deliver 
aX P tan vr ( n ) as the generated value. 


3.15 Irwin-Hall distribution 

Let Uk ^ W(0, 1) fc = l,...,n n uniform in- 

n 

dependent distributions. In this case, X = '^^Uk 

k=l 

follows the Irwin-Hall distribution with parameter n 

{X-^in{n) ). 

Its probability density function is given by: 


fx{x) = 


1 


2(n- 1)! 


E(-i) 


k=0 


k (n 

k 


(x — ky ^ sgn(x — k) 



with X G [0, n] and sgn(x) the function given by: 

— 1 if X < 0 

sgn(x) = ■( 0 if X = 0 

1 if X > 0 


0.43636 0.87273 1.3091 1.7455 2.1818 2.6182 3.0545 3.4909 3.9273 4.3636 4.8 


in(5) 
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To generate X ^ 


return 


n 

x = 'Y^Uk 

k=l 


IH{n), Ui,U2,...,Un 
as the generated value. 


values must be generated from U{0,1) 


and then 


4 Discrete distributions random generation 


In this section different procedures are presented in order to generate discrete distributions. The 
inverse transform method for discrete distributions described in section 2T is used in order to generate 
a sample of the distribution. The only thing to do is to use the random_discrete function developed 
in the same section with the corresponding parameters. 


4.1 Uniform discrete distribution 

X = {o, a+1 , . . . ,b} has an uniform discrete distribution with parameters a and b { X ^ b) ) 

if its distribution mass function F is: 

Fix) = P\X = x] = ; X = a, a + 1, . . . , 6 

^ ^ b-a+1 

Thus, in order to generate a sample of size n from an Uv{a, b), the following Derive code can 
be used: 


random_unif orm_discrete (n :=1, a:=0, b :=1) : = 
random_discrete(n, 1/ (b-a+1) ,a) 


4.2 Bernouille distribution 

X = {0, 1} has a Bernouille distribution with parameter p { X ^ Bev{p) ) if its distribution mass 
function F is: 

F{x) ^ P[X ^ x]=p^{l-p)^-^ ; x = 0,l 

or, equivalently, F{0) = 1 — p and T(l) = p. 

Thus, in order to generate a sample of size n from a Bev{p), the following Derive code can 
be used: 

random_bernouille(n := 1, p := 1/2) := random_discrete(n,p^(l — p)^~^,0) 

4.3 Rademacher distribution 

X = { — 1,1} has a Rademacher distribution ( X ^ TZad ) if its distribution mass function F is: 

F{x) = P[X = x\ = ^ ; X = —1, 1 

If T Ber(l/2) then X = 2y - 1 TZad. Thus, in order to generate a sample of size n 
from a TZad, the following Derive code can be used: 

random_rademacher (n := 1) := VECT0R(2k-l ,k,random_bernouille(n, 1/2) ) 
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4.4 Binomial distribution 


X = has a binomial distribution with parameters n and p ( X Bi{n,p) ) ii its 

distribution mass function F is: 

F(x) = P[X = x]= ^ ^p^{l-pr-^ ; x = 0,l,...,n 

Thus, in order to generate a sample of size m from a Bi{n,p), the following Derive code can 
be used: 


random_binomial(m := 1, n := 100, p := 1/2) := 
random_discrete(m, ncom(n,x) p^(l — p)“~’',0) 


4.5 Poisson distribution 

X = (0, 1,2,.. .} has a poisson distribution of parameter A { X ) if ifs distribution mass 

function T is: ^ 

F{x) = P[X = x] = ^ ^ ; x = 0,l,2,... 

X . 

Thus, in order to generate a sample of size n from a P(A), the following Derive code can be 
used: 


random_poisson(n := 1, A := 1) 


random_discrete (n 



, 0 ) 


4.6 Geometric distribution 

X = {0,1,2,...} has a Geometric distribution with parameter p ( T ^e(p) ) if its distribution 
mass function F is: 

F(x) = P[X — x] = p {1 — pY ; X = 0, 1, 2, . . . 


Thus, in order to generate a sample of size n from a Ge{p), the following Derive code can be 
used: 


random_geometric(n := 1, p := 1/2) := random_discrete (n,p (1 — p)^,0) 


4.7 Negative Binomial distribution 

X = {0, 1 , 2 ,.. .} has a negative binomial distribution with parameters r and p { X ^ J\fB{r,p) ) 
if its distribution mass function F is: 

F{x) = P[X = x]= ^^pYI-pY ; a; = 0,1,2,... 

Thus, in order to generate a sample of size n from a J\fB{r,p), the following Derive code can 
be used: 


random_negative_binomial(n := 1, r := 10, p := 1/2) : = 
random_discrete(n, ncom(r+x-l ,x) p’^(l — p)’',0) 
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4.8 Hypergeometric distribution 

X = (max(0,ni + m — n), . . . , min(n,i, m)} has an hypergeometric distribution with parameters 
n , m and ni {X ^ 7i{n,m,ni) ) if its distribution mass function F is: 


F{x) = P[X = x = 


rii \ I n — rii 
X } \ m — X 

n 
m 


Thus, in order to generate a sample of size 
be used: 


X = max(0, rii + m — n), , min(ni , m) 


from a J\fB{r,p), the following Derive code can 


random_hypergeometric(l:=l, n := 100, m := 50, nl := 50) := 

ncominl, x)ncom(n — nl, m — x) 

random_discrete(l, 7 r ,max(0, nl + m - n)) 

ncom(n, m) 


5 Approximative algorithms 

In this section several algorithms used to generate different distributions by approximations are 
presented. 

The main reason in order to use such algorithm is that they produce “good” samples and they 
are quite much faster than other “exact” algorithm. In fact, some of the algorithms described 
above are approximative algorithm (random_gainma9, random_beta7 and random_chi_square2 are 
approximative algorithms for Be(a, P) and distributions which use 

distribution as approximation). 

In the followings subsections some other approximative algorithms which have been implemented 
in Random_distribution package are developed. 

5.1 Approximation to Binomial distribution by Poisson distribution 

Let X ^ Bi{n,p). If p is “small” then 


Bi{n,p) ~ F{np) 

Thus, the algorithm to generate a sample of size m from the binomial distribution approximated 
by a poisson distribution is: 


random_binomial_approx_poisson(m:=l, n := 100, p:= 0.01) : = 
random_poisson (m , np) 


5.2 Approximation to Binomial distribution by Normal distribution 

Let X Bi{n,p). If n is “large” and p is not close to 0 or 1, then 

Bi{n,p) Ff (^np, \/np (1 — p)^ 

It can be shown that “good” approximations are reached when n p > 10 with 0 << p < 0.5 or 
n (1 — p) > 10 with 0.5 < p << 1 ( << indicates less than but not close to). 

On the other hand, as Bi{n,p) has only nonnegative integer values, the nearest nonnegative 
integer to the value returned by Ff (np, y/np (1 — p)') must be chosen as the generated value for 
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X. It is easy to understand that the nearest nonnegative integer to a value z can be found by the 
operation: 

max ^0, floor 

where floor(x) is the integer part of x. 

Thus, the algorithm to generate a sample of size m from the binomial distribution approximated 
by a normal distribution is: 



random_binomial_approx_normal(m:=l, n := 100, p:= 1/4) : = 

vector (max (0 , floor (k+1/2) ) , k , randommormal (m , np , sqrt (np ( 1-p) ) ) ) 


5.3 Approximation to Poisson distribution by Normal distribution 

Let X^V{\). If A >10 then V{\)^M 

On the other hand, as P(A) has only nonnegative integer values, the nearest nonnegative integer 
to the value returned by A/" ^A, must be chosen as the generated value for X. As shown in 


section 5.2, this can be easily done by the operation: 


max ^0, floor + 2 

Thus, the algorithm to generate a sample of size n from the poisson distribution approximated 
by a normal distribution is: 


random_poisson_approx_normal (n := 1, A:= 15) : = 

vector (max(0 .floor (k+1/2) ) ,k , random_normal (n, A , sqrt (A) ) ) 


5.4 Approximation to Geometric distribution by Exponential distribu- 
tion 


li X Qe{p) then X ^ \ In 


1 


-P 

On the other hand, as Qe(p) has only nonnegative integer values, the nearest nonnegative 
integer to the value returned by £ fin f — - — ] ) must be chosen as the generated value for X. 


1 — p j 

The nearest nonnegative integer to a value z > 0 (exponential distribution has only positive real 
numbers) can be found by the operation: 

floor ( ^ + 2 


where floor(x) is the integer part of x. 

Thus, the algorithm to generate a sample of size n for the Geometric distribution approximated 
by a exponential distribution is: 


random_geometric_approx_exponential(n := 1, p:= 1/2) : = 

vector (floor (k+1/2) , k , random_exponent ial (n,ln(l/(l-p)))) 
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6 Graphical Approach 

In order to check graphically if the generated sample fix the distribution, two different graphical 
functions have been developed. 

• bars_diagram(v, k, mi, ma): it returns the needed instructions to plot the bars diagram 

of vector v divided in k classes (k=20, by default). The optional parameters mi (minimum) 
and ma (maximum) are the minimum and maximum values of vector v unless other values are 
considered. This values will be used to set the plot range. 

• density_f unction (name, para): it returns the corresponding density function of distri- 

bution name which parameters are set in vector para. The possible values for name are: 
“Uniform”, “Exponential”, “Normal”, “Lognormal”, “Cauchy”, “Weibul”, “Gamma”, “Beta”, 
“Chi-square”, “t”, “F”, “Z”, “Pareto”, “Logistic” or “IrwinHall”. 

Note that these functions have been developed for continuous distributions. 

7 Examples 

As an example, let us describe the steps needed to obtain a sample of size n = 1000 from a Normal 
distribution with parameters p = 3 and a = 2: 

1. Load the utility file Random_distributions .mth in DERIVE 6.1. 

2. Write down the following instruction in the Author Line: random_normal(1000,3,2) 

3. Do not simplify but Approximate the expression (in order case, time would be too long). 

4. Variable res stores the result. 

Now, in order to check graphically if the generated sample fix the distribution, the steps to 
follow are: 

5. Approximate the expression: bars_diagram(res) 

6. Follows the instructions given after this execution (the appropriate range for the plot) and plot 
the result. 

7. Simplify the expression: density_funciton("Normal" , [3,2] ) and plot the result. 

Remind that res variable stores the sample (of size n — 1000 in the example). If, the sample 
does not fix the density function, a new sample can be generated again going to step 1. 

In order to see more examples of different generated samples from discrete and continuous dis- 
tributions, as well as, examples of different plots, please refer to files Random_distributions .dfw 
and Random_distributions .mth. These files are distributed together with this paper and are also 
allocate at http : / /www . matap . uma . es/prof esor/ j l_galan/Derive/Packages. 


8 Future work 

In this section some possible future work lines are described: 
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• Development of different programs to generate other discrete or continuous distributions: al- 
though in this work the most common distribution has been considered, there exist many more 
distributions which can be generated. This work shows different ways of generating samples 
for both discrete and continuous distributions. Even more, since some generic algorithm has 
been developed they may be used to generate almost any new distribution (mainly for discrete 
distributions) . 

• Development of graphical approaches for discrete distributions: since the graphical approach de- 
veloped has been optimized for continuous distributions, some graphical approaches for discrete 
distributions could be also considered in future. 

• Development of applications: there exist many applications which requires working with sam- 
ples of different distributions. This kind of applications can be considered using this work. 

• Checking the obtained result theoretically: although a graphical approach has been developed to 
check if the generated sample fix the corresponding distribution, different theoretical approach 
can be considered in order to check the results. 

• Translation of these algorithms: Although this package has been developed in Derive, it can 
be easily translated to other Computer Algebra Systems (Cas) or even to a portable platform 
as Java. 
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Derive- and CAS-TI-User Forum 


D-N-L#75 


From: DERIVE computer algebra system [mailto:DERIVE-NEWS0JISCMAIL.AC.UK] On 
Behalf Of Volker Loose 

Sent: Thursday, October 08, 2009 9:10 AM 
To : DERIVE-NEWS0 JISCMAIL.AC.uk 
S ubject: nsolve 


nsolve gives only 1 solution of an equation. I tried to write a function 
vnsolve that gives more solutions, using the vector- and the nsolve- func- 
tion. It seems to be a mistake in my function but I can't find it. Has anyo- 
ne an idea? 


Many thanks 
Volker 

X 

#1: f(x) \- e. - X - 2 


#2: vnso1ve(f, u, o, n) VECTOR 


- 1) ■ (o - u) 


1 ■ (o - u) 


, 1 , 1, n, 1 


(i - !)■ (o - u) I- (o - u) 

NSOLVEl f-x, X, u + , u + 


#3: vn5o1ve(f (x) , -B , B, 2) 



1 


2 



-B 


0 



0 

1 

B 



X, 

II 

o 


. X = 0 . 



DNL: Do you find the bug? Volker defined in #1 function f(x) and then in #2 he used f as a parameter 
for his function. Please note in the 4"^ component of the vector nsolve ( f x, ...) . This should be f(x) 
but is interpreted as a product f x. Replace f in the parameter list and f x in the NSOLVE-command by, 
say f_, then vnsolve works. 

But if there are more solutions which are close then this procedure does not work properly. 

Try vnsolve (SIN (x^2) -e^x, -7,7,10) which should give 16 zeros, vnsolve gives only 5. 

Jim FitzSimmons suggested another procedure: 





r fCx) 7 



ITERATE 

APPEND 

NSOLUTIONS 

. X. -B, B 

^ n(x - y, y, v) ^ 

1 V 

J 

, V, [], 2 

J 


[1.146153220, -1.841405660] 


In case of more and close zeros it does not a satisfying job, too. In DNL#63 I tried to define a function 
finding the zeros. See mysolutions: 

2 X 

mysoluti'onr2(SIN(x ) - e , x, -7, 7) 

[-6.86460S15S, -6.632014374, -6.390542007, -6.140135705, -5 .S7S32626S , -5.605315329, -5.316900028, -5.013915234, -4.688491003, 
-4.343104006, -3.960923705, -3.548961687, -3.062349296, -2.522602266, -1.720968215, -0.7149689691] 


At the occasion of this exchange of e-mails Danny Ross Lunsford wrote: 
about DERIVE-NEWS(3>JISCMAIL.AC.UK 


I am very happy to see that this list still exists. Derive is still a nearly ideal math tool and I 
hope it lives on. 

I have nothing to add, except this: I completely agree with Danny, Josef. 


